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SUMMARY 

The capability of the boundary element method (BEM) in determining thermal boundary 
conditions on surfaces of a conducting solid where such quantities are unknown has been 
demonstrated. 1 he method uses a non-iterative direct approach in solving what is usually called the 
inverse heat conduction problem (IHCP). Given any over-specified thermal boundary conditions 
such as a combination of temperature and heat flux on a surface where such data is readily available, 
the algorithm computes the temperature field within the object and any unknown thermal boundary 
conditions on surfaces where thermal boundary' values are unavailable. A two-dimensional, stead v- 
state BEM program has been developed and was tested on several simple geometries where the 
analytic solution was known. Results obtained with the BEM w'ere in excellent acreement with the 
analytic values. The algorithm is highly flexible in treating complex geometries, mixed thermal 
boundary conditions and temperature-dependent material properties and is presently being extended 
to three-dimensional and unsteady heat conduction problems. The accuracy and reliability of this 
technique was very good but tended to deteriorate when the known surface conditions w’crc oniv 
slightly ovcr-spccified and far from the inaccessible surface. 


INTRODUCTION 

The objective of the steady-state inverse heat conduction problem is to deduce temperatures 
and heat fluxes on any surface or surface element where such information is unknown. In manv 
instances it is impossible to place sensors and take measurements on a particular surface of a 
conducting solid due to the inaccessibility or seventy of the environment on that surface These 
unknown thermal boundary values may be deduced from additional temperature or heat flux 
measurements made writhin the solid or on some other surface of the solid. This problem has been 
given a considerable amount of attention by a variety of researchers and virtually all w'ork has been 
directed to the one-dimensional transient problem. The first method proposed to solve the IHCP 
used inversion of convolution integrals (Stoiz 1960) and was subsequently improved bv a number 
of authors (Beck et al. 1988). Many other methods have also been developed using such 
techniques as Uplace transforms, finite elements, time-marching finite differences and other 
099?) CheS ‘ A detaiied chronol °g lcal review of the IHCP literature has been provided by Hensel 

... . A characteristic of most of these inverse techniques is that they tend to produce temporal 
oscillations in the unknown surface thermal condition estimates that are larger than the temporal 
osciUauons in the over-specified thermal data as it propagates through the solid (Hills and Hensel 
1986). In other words, the random noise due to round off errors tends to magnify as the solution 
proceeds and quickly produces a useless solution, especially as the distance between the surface 
and the over-specified information increases. A number of authors have presented v arious 
smoothing techniques for reducing this error growth, but the effect of these operations on the 
accuracy of the solution is not easy to evaluate (Murio 1993). 
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The method presented herein doesrnot utilize any artificial smoothing technique and is not 
limited to transient or one-dimensional problems. This approach is non-iterative and has been 
shown to compute meaningful and accurate thermal fields in a single analysis using a straight- 
forward modification to the boundary element method (BEM). 6 

The BEM is a very accurate and efficient technique that can solve boundary value problems 
such as those governing heat conduction, electromagnetic fields, irrotationai incompressible fluid 
flow, elasticity and many other physical phenomenon. For steady-state heat conduction analysis 
using the BEM, either temperatures, T, or heat fluxes, Q, are specified everywhere on the surface of 
the wild where one of these quantities is known while the other is unknown. In the BEM solution 
to the IHCP, both T and Q must be specified on a part of the solid's surface, while both T and O are 
unknown on another part of the surface. Elsewhere on the solid's surface, normal boundary 
conditions should be applied as either T's or Q's. The surface section where both T and O 'are 
specified simultaneously is called the over-specified boundary and is necessary for the IHCP 
problem's solution. } 

nmhi^rr^* ca ! multipl y connected, inverse heat conduction 

problem. Surfaces labeled T\ are the over-specified boundaries where both T and Q are civen 

Normal boundary conditions (either T or Q specified) are enforced on the surfaces labeled Ti. 
Thermal data is assumed to be inaccessible on the inner T 3 surface and thus has both T and 6 

on this boundary. The objective of the IHCP is to compute temperatures and heat fluxes 
on the boundary' T 3 using only the values of T and Q provided on the surface of the solid and 
possibly, additional temperature measurements made within the solid if such data is available ’ 


THEORY 

Two-Dimensional Steady-State BEM 

,, Steady-state heat conduction in a homogeneous medium with a constant coefficient of 
thermal conductivity is governed by the Laplace's equation in the region, Q, of a conducting solid 

v2 T = 0 (1) 

n'lfZlt ? e ter "P eratur fJ hi s is a linear boundary value problem having essential boundary' 
conditions, T 0 , and natural boundary condiuons, Q 0 , specified on the surfaces Tu and Tq > 

^Sntgive 0 n r b n y 0nlmear ** “^^P^ent matenal properties, the governing 


V* (HT) VT) = 0 


(2) 


where X(T) is the temperature-dependent thermal conductivity. Equation 
the application of the classical Kirchoff transformation which defines the 


(2) can be linearized 
heat function, ©, as 


by 


e = j HE! dT 


0 ^0 


(3) 


Here, ^ is a reference conductivity and X(T) could be an arbitrary function of temperature. 
Consequently, equation (2) can be transformed into Laplace's equation and solved for the heat 
£nct,on, ©, instead of temperature, T. Results obtained for the heat function must £ transited 
back into temperatures using the inverse of the transformation given in equation (3). 
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Laplace's equation may be solvednising the BEM (a weighted residual technique) by 
introducing an approximation, u, to the exact solution, 9. Since the approximation is, in general, not 
equal to the exact solution, an error function or residual is produced in the domain and on the 

boundary. The residual in the domain is given by R = V 2 u and the residuals at the boundaries are 
R u = u - ©o and R q = du/dn - Qo< These error functions are normally non-zero unless u is the 
exact solution. The weighted average of the residual over the domain and on the boundary may be 
set to zero by the weighted residual statement 


/ u* V 2 u dQ - f (q - Q 0 )u* dT + / (u - 0 o )q* dT = 0 (4 

a r q r\i 

where u* represents the weight function which is usually called the fundamental solution (Brebbia 
and Dominguez, 1989), while q = du/<3n, q* = du*/dn and n is the direction of the outward normal 

to the surface T. After integrating by parts twice, the boundary integral equation for Laplace's 

equation is obtained 


f u V- u * d£2 + f u* q dT = f q * u dT 

a t r 

The weight function is a Green's function solution for a point-source subject to the 
homogeneous boundary conditions. For the two-dimensional Laplace's equation it is 


"* = 2jt'° 7 8 (r) <«> 

where r — I x; - Xj |, xj is the coordinate of the observation point, xj is the coordinate of the source 
point and the logarithm function here has base e. The bounding surface V is discretized into N 
surface elements bounded by N end-nodes. After discretizing the surface and utilizing the 
properties of the Dirac delta function, the boundary integral equation (6) can be written as 


CjUJ + 



dTj 


(7) 


for each ith node. The term q indicates the scaled internal angle at the ith surface node, 
junctions u and q are assumed to vary linearly along each surface element and, therefore 
be denned in terms of their nodal values and interpolation functions 


The 

they can 


u (!) = 4>1 (!) u i + fc(D u 2 and q (|) = <|>i(|) q i + ^d) q 2 ( 

where | is a localized surface-following dimensionless coordinate, while <j>j - (l - £)/2 and 

“ (1 + !)/^* T'he whole set of equations for the N nodal values of u and q can be expressed in 
matrix form as K 
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m 'U * [G] Q 


(9) 


where U (Ui,U 2 ,U 3 ,...,Un) und Q — (Qi>Q2»Q3»" 4 »Qn) ure vectors containing the nodal potentials 
and surface panel fluxes while the terms in the [H] and [G] matrices are assembled bv properly 
adding the contributions from each surface integral 


Hjj = f<M* <Hj + J>, q* dr^j 

r >, 


G|J = J> 2 u* drj + J>,u* dT j+1 
r j 

The free term, q, is produced when the first surface integral of equation (6) is 
integrated in the sense of the Cauchy principal value. Since q* = du*/6m = (3u*/dr) (dr/dn) = 0 
when the ith surface integral contains the ith observation point, the diagonal of the [H] matrix is 
amply the q term. This coefficient may be computed explicitly by calculating the internal angle at 
the surface node or implicitly (Brebbia and Dominguez 1989) by first assuming a constant unit 
potential throughout the enure domain and then solving for the diagonal component as 


N 

c i = H ii = - 2 H ij ' * J (11) 

When the observation node is on the surface panel of integration, the terms in the TGI matrix are 
computed analytically from the integral 


G ” ~ 2^/*2 i0 §(r) dr ' + 2^ Ai ,0 s(~) dr i 


+1 


(12) 


i+l 


After the [H] and [G] matrices are formed, all boundary conditions arc applied and a set of 
linear algebraic equations [A] X = F, is constructed. Known or specified surface potentials, U, 
and fluxes, Q, are assembled on the nght-hand-side of the equation set and are multiplied by tiieir 
respective [H] or [G] matrix row thus forming the vector of knowns, F. All unknown potentials or 
fluxes are assembled on the left-hand-side of the equation set and are represented by a coefficient 
matrix [A] multiplying a vector of unknown quantities, X. 

_ The i et of linear algebraic equations is then solved for the unknown surface potentials U 
and fluxw Q, using a singular value decomposition (SVD) matrix solver (Press et al 199^) We 
used the SVD since the matrix tends to become ill-conditioned or singular (several equations 
become linearly dependent with other equations in the equation set) whenever the over-specified 
thermal data are farther away from the surfaces where no boundary conditions are applied If 
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additional thermal data within the solid irprovided, additional equations may be added to the 
equation set. Note that the SVD algorithm is capable of providing a satisfactory solution vector 
even when the [A] matrix is not square. The more rows (i.e. more data points) that are provided to 
the system, the more accurate the solution vector becomes, although the reverse is true when the 
matrix has less rows than columns. Once the matrix is solved, the entire thermal field within the 
solid can be easily deduced. 


RESULTS AND DISCUSSION 
IHCP for a Square Plate Using the BEM 

A BEM computer program was developed using the theory discussed in the previous 
section. The accuracy of the BEM as a solution to the IHCP was verified for a solid square plate. 
The plate was 6.0 m on each side and the thermal conductivity of the plate was chosen as 1.0 W/mK. 
The top and bottom boundaries were specified to be adiabatic (Qo = 0 W/m 2 ) while the left side 
surface of the plate was over-specified with a temperature boundary condition of To = 300 K and a 
heat flux boundary condition of Qo = -50 W/m 2 . The right side boundary was considered to be 
inaccessible and, as such, both temperature and heat flux were unknown on this boundary. The 
plate boundary was discretized with 12 panels (3 per each of the four sides) on the boundary of the 
solid. The BEM was successful in computing a temperature field within the plate that was accurate 
to almost the floating point precision of the computer. The computed temperature and heat flux on 
the right side boundary were 0.000161 K and 49.99997 W/m 2 , respectively. 

Study of the IHCP for an Annular Disk Using the BEM 

The behavior of this algorithm for various combinations of boundary conditions was 
documented for steady-state heat conduction in an annular solid disk. The outer radius of the disk 
was 1.2 in and the centrally located hole had a radius of 0.5 in. The analytic solution for this 
problem was developed by applying Dirichlet or essential boundary conditions everywhere on the 
boundary of the annular region. Temperature boundary conditions of 100°C on the outer boundary 
and 50°C on the inner boundary were enforced. The thermal conductivity of the solid was 
considered to be constant, X — 1.0 Btu/in sec°R. The analytic solution for the temperature field 
within the disk is easily found as 


T(r) = A + Blogr ( 13 ) 

where A = 89.59 and B = 57. 1 1. The radial heat flux is then 

Q (r) = -X VT = - X dT(r)/dr = B / r ( 14 ) 

which yields Q out = -47.59 Btu/in 2 sec and Qi n = 1 14.22 Btu/in^ec as heat fluxes through the 
outer and inner boundaries, respectively. The BEM algorithm was run on the same problem The 
problem was discretized with 36 panels on each outer and inner boundary'. The BEM program 
predicted the temperature field in the solid which averaged only a 0.3% error versus the analytic 
solution. 

In order to study the feasibility and accuracy of the BEM solution to the steadv-state IHCP 
seven variations to the same problem were performed and the results obtained were compared to 
those from the previous problem. Each test utilized the same annular geometrv and outer boundary 
thermal data in a variety of combinations. 

Tm L The outer and inner boundaries of the annular domain were each discretized with 36 
equally-length fiat panels. The entire outer boundary was over-specified with temperature and flux 
boundary conditions, while both temperature and flux were unknown on the inner boundary'. The 
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BEM formulation detailed in the theory action, can be represented in matrix form by equation (9). 
For this test case, the solution set of 72 equations included 72 known values given as boundary 
conditions on the outer surface and 72 unknowns on the inner boundary'. The BEM computed the 
temperature Held within the annular solid in addition to the unknown temperatures and heat fluxes 
on the inner boundary. Figure 2a shows the computed temperature contours for the annular solid 
disk and also includes the BEM nodes used and the type of boundary conditions specified at each 
node. The box-shaped nodes have both T and Q known and thus are* over-specified, the circles are 
nodes where both T and Q are unknown, and the triangular nodes have a single boundary condition 
of temperature applied. The thick solid lines in figures 3 and 4 represent the accuracy of this 
particular BEM solution. Figure 3 shows the relative percentage error in temperature on the inner 
boundary for each test as a function of the circumferential angle in radians. Figure 4 is the same as 
figure 3 except that it gives the relative percentage error in the heat flux on the inner boundary. 
Notice that this test case had an almost perfectly symmetric result with an average error of only 
0.5% in temperature and a somewhat oscillating error in heat flux averaging about -1.5%. 

Test 2. This test case was identical to Test 1 except that the outer and inner boundaries were 
discretized with a coarser grid consisting of 18 panels each. Overall, the BEM solution set had 36 
knowns, 36 unknowns and 36 equations. The computed temperature field and boundary 
discretization are shown in Figure 2b and the relative percentage error in temperature and heat flux 
on the inner boundary are given as thin solid lines in figures 3 and 4. The temperature field within 
the solid was nearly perfectly symmetric, but was uniformly biased about 2.5% in temperature and - 
2.5% in heat flux. 

Test 3. This test case was identical to Test 1 except that the boundary of the annular disk 
was discretized with 36 panels on the outer boundary' and 18 panels on the inner boundary. 

Overall, the BEM solution set was over-specified and had 72 knowns, 36 unknowns and 54 
equations. The thick dotted lines in figures 3 and 4 readily show that Test 3 produced the most 
accurate results for both temperature and heat flux. In addition, the temperature contours in Figure 
2c are nearly perfectly symmetric. 

T e$t 4. This test was identical to Test 1 except that the outer boundary' was discretized with 
18 panels and the inner boundary was discretized with 36 panels. Overall, the BEM solution set 
was under-specified and had 36 knowns, 72 unknowns and 54 equations. Results of Test 4 are 
given by the thin dotted line in figures 3 and 4. The temperature was uniformly biased with a 3.0% 
error, while the heat flux was somewhat oscillatory and similarly biased. The temperature contours 
in figure 2d were nearly symmetric. 

Test 5. Both the outer and inner boundaries of the annular disk were discretized with 36 
panels. Temperature boundary conditions were specified everywhere on the outer boundary but the 
additional heat flux boundary conditions were over-specified in the first and third quadrants of the 
outer boundary only. The BEM solution set had 54 knowns, 90 unknowns and 72 equations The 
temperature field shown in figure 2e was comparable to that of Test 4. The temperature and heat 
flux on the inner boundary are represented by the finely dotted lines. The temperature distribution 
on the inner boundary was somewhat oscillatory, but averaged only a 0.75% error. The heat flux 

on the inner boundary was also oscillatory' and averaged an error of about - 2 . 0 %. 

Is st 6. The circular disk was discretized with 36 panels on both the*inner and outer 
boundaries. Temperature boundary conditions were specified on the entire outer boundary while 
heat flux boundary conditions were over-specified only on the upper half of the outer boundary. 

As in Test 5, the BEM solution set contained 54 unknowns, 90 unknowns and 72 equations The 
temperature field illustrated in figure 2f was asymmetric about the x-axis, but was very nearly 
symmetric about the y-axis. The greatest error in the temperature field occurred in the bottom half 
of the annular solid region. The thick dashed lines in figures 3 and 4 reveal inner boundary errors 
that are quite oscillatory in nature and noticeably peak at the very bottom of the solid disk (the point 
farthest from the over-specified data). ^ 

T est li This test case is identical to Test 5, except that heat flux boundary conditions are 
over-specified in the first quadrant of the outer boundary only. The BEM solution set contained 45 
UI ^ cnowns ^ T2 equations. Figure 2g illustrates the temperature contours within the 
solid disk. The error in the temperature field obviously worsens as the distance from the over- 
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specified data increases. The thin dashecflines in Figures 3 and 4 illustrate the error in the 
temperature and heat flux on the inner boundary. The error is oscillatory and peaks at about 60% at 
the point farthest from the over-specified data. Notice also that the temperature field is symmetric 
about the line inclined 45 degrees and passing through the center of the circle. 

IHCP for a Rocket Nozzle Wall Section with a Coolant Flow Passage 

The BEM solution to the IHCP was attempted on a realistic engineering problem with 
temperature-dependent material properties. High pressure, reusable rocket thrust chambers 
encounter a progressive thinning of the coolant flow passage wall after repetitive engine operation. 
This deformation is caused by high thermal plastic strains that eventually cause cracks to form in 
the cooling passage wall. An engineer who wishes to reduce or eliminare the plastic strain may 
obtain experimental data such as hot gas wall temperatures and heat fluxes, shroud temperatures, 
compressive strains, and thrust chamber total pressure and temperature (Quentmeyer 1978, 1992). 
Unfortunately, the engineer cannot normally obtain data within the coolant flow passage due to the 
extremely low temperature of the liquid hydrogen coolant and the small dimensions of the passage.. 

Figure 5 is a schematic of a cylindrical thrust chamber assembly and figure 6 illustrates a 
cylinder wall cross section showing typical instrumentation locations and dimensions. These 
figures were taken from a NASA publication (Quentmeyer 1978) and were subsequently used to 
generate the geometry of the nozzle wall section. The hot gas wall temperature (1520 °R), heal flux 
(-35 Btu/in 2 sec) and shroud temperature (518.4 °R) were experimental measurements taken from 
the same publication. The outer shroud heat flux was assumed to be negligible (0 Btu/in 2 sec). The 
coefficient of thermal conductivity of the solid copper region was linearly dependent on the local 
temperature 

X. - Xq ( 1 + a T ) (15) 

where Xo - 0.004893 Btu/insec°R and a = -0.000055056 °R-*. In addition, figure 6 shows that the 
conducting solid region is made up of three different materials; copper, electrodeposited copper, and 
nichrome ZrC> 2 . Although the present analysis uses only a single material, the BEM can be 
modified to handle composite materials with each having distinct thermal properties. The 

shaded portion in figure 6 is the domain typically used in the two-dimensional heat conduction 
model. For the BEM analysis of this IHCP, a full section containing the entire cooling passage 
and half of the surrounding conducting metal was generated in order to examine the symmetry' of 
the results. The meridional or symmetry planes were assumed to be adiabatic. The outer andinner 
boundaries were discretized in the same manner 16 panels on the hot gas side , 8 panels on the 
shroud and 8 panels on each of the two periodic meridional boundaries. The BEM solution set 
contained 66 knowns, 94 unknowns and 80 equations . The BEM computed both temperatures and 
heat fluxes on the entire coolant flow passage boundary in addition to the temperatures on the 
meridional side boundaries. The predicted temperature field within the solid region is illustrated in 
figure 7. These results show a negligible asymmetry about the meridional centerline and slight 
oscillations in the temperatures computed near the comers of the coolant passage's cool side. 


CONCLUSIONS 

The boundary element method computed temperature and heat flux boundary conditions on 
boundaries of a conducting soiid where such quantities were originally inaccessible and unknown 
The results presented herein indicate that the direct non-iterative BEM* solution method for the 
IHCP is an accurate, robust and reliable technique that takes only seconds of CPU time on any 
typical mainframe, workstation or PC. In addition, the results obtained were found to be more" 
accurate when one or both of the following conditions were observed: a) greater amount of over- 
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specified data was applied, b) the over-specified data locations were in close geometric proximity to 

the locations of the unknown boundary conditions. F y 


ACKNOWLEDGEMENTS 

Sincere thanks are due Mr. Vineet Ahuja for suggesting the use of the SVD algorithm for 
almost singular matrices since a standard matrix solver produced meaningless results due to severe 
accumulation of round-off errors. 


REFERENCES 


Beck ’ B. and St Clair, CR., Jr.: Inverse Heat Conduction: Ill-Posed Problems 

Wiley-Interscience, New York, 1985. 

Brebbia, C.A. and Dominguez, J.: Boundary Elements, An Introductory Course. McGraw- 
Hill Book Company, New York, 1989. 

Hensel, E.C., Jr.: Multi -dimensional Inverse Heat Conduction. Ph.D. dissertation 

Mechanical Engineering Dept., New Mexico State University, Las Cruces NM 1986 

Hills, RjG. and Hensel, EC., Jr.: One-dimensional Nonlinear Inverse Heat Conduction 
Technique. Numerical Heat Transfer, vol. 10, pp. 369-393, 1986 

Martin, T.J.: Inverse Design and Optimization of Two and Three-Dimensional Coolant Row 
Suvosfty ^ay T 1993 S ’ DepL ° f Aeros P ace Engineering, The Pennsylvania State 

Murio, D.A.: The Mollification Method and the Numerical Solution of Ill-Posed Problems 
John Wiley & Sons, In., New York, 1993. 

Press, W^ Teukolsky, S.A. Vetteriing, W.T. and Hannery, B.P.: Numerical Recipes in 
FORTRAN. Second Edition, Cambridge University Press 199'’ 

NASA TM-^S ^ ° f °" R ° Cke ' Thnlst Chamber 

Quenone^RJ^ An^pen^nu 1 Investigation of High-Aspect-Ratio Cooling Passages. 

St olz, G. Jr.: Numerical Sol ulions lo an Inverse Problem of Heal Conduction for Simple 
Shapes. ASME Journal of Heat Transfer, vol. 82, pp. 20-26, 1960. 



Figure 1. A geometric definition of a two-dimensional inverse heat conduction problem. 
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Geometry of the BEM nodes on the outer and inner boundaries, boundary 
condiuon types and isotherms computed with the BEM for each of the sev 
annular disk test cases. 
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Figure 3. Relative percentage errors (BEM versus analytic solution) of the inner 
boundary temperatures for each of the seven annular disk test cases. 
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Figure 4. 


Relative percentage errors (BEM versus analytic solution) of the inner 
boundary heat fluxes for each of the seven annular disk test cases 
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1 520 °R, Q * -35 Btu / in 2 see) 



(T Unknown, Q = 0 Btu / in 2 see) 


Figure 7. 
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518 4 °R. Q = 0 Btu / in 2 see) 



